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ABSTRACT 


The  time-averagad  power-density  distribution  below 
the  ocean  surface  due  to  incident  laser  radiation  is 
examined  by  means  of  computer  simulation  of  the 
geometrical  optics  involved  with  the  air/sea  interface 
and  subsequent  ocean  penetration  by  the  laser  beam. 
The  effects  over  the  entire  spectra  cf  incidence 
angles,  wind  velocities,  wind  directions,  beam  spot 
sizes  and  depths  of  penetration  are  analyzed. 
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I.   INTRODUCTION 


The  time-averaged  power-density  distribution  below  the 
ocean  surface  due  to  incident  laser  radiation  is  examined 
over  the  entire  spectra  of  incidence  angles,  wind  velocities 
and  wind  directions  by  means  of  geometrical  optics. 

Current  interest  in  using  laser  systems  for 
communications  through  the  air-sea  interface  and  for 
detecting  submerged  objects  necessitates  the  development  of 
prediction  aethods  for  determining  the  subsurface 
power-distribution  of  radiation  from  a  laser  source  above 
the  ocean.  A  mathematical  model  for  the 
optical-communications  application  done  by  Karp  in  Hef.  1  is 
developed  in  terms  of  a  radiance  function  related  to  the 
mutual  coherence  function  with  accurate  results  for 
incidence  angles  out  to  45  degrees.  A  more  general  model 
developed  by  Swennen  in  Refs.  2-4  examines  the  power-density 
distribution  below  the  ocean  surface  through  a  rigorous 
assessment  cf  the  ocean  surface  geometry  allowing  for 
theoretically  accurate  results  over  the  entire  incidence 
angle  spectrum. 

The  computer  code  developed  for  this  analysis  of  the 
power-density  distribution  below  the  ocean  surface  is  an 
expansion  of  Swennen's  theory[Refs.  2-4]. 


II.   DESCRIPTION  OF  PHYSICAL  PARAMETERS 


A.   OCEAN  SURFACE 


A  geometrical  representation  of  the  ocean  surface  in 
terms  of  local  slopes  developed  by  Cox  and  Munk[Ref.  5]  is 
the  basis  for  the  analysis  of  the  power  probability 
distribution  telow  the  ocean,  as  it  is  in  Refs.  2-4.  A 
review  of  this  representation  follows. 

The  center  of  symmetry  of  the  incident  laser  beam  en  tne 
ocean  surface  facet  is  the  center  of  a  right-  handed 
cartesian  coordinate  system,  point  "0"  in  Fig.  1.  The 
z-axis  points  vertically  upward.  The  y-axis  is  in  the 
horizontal  plane,  colinear  with  the  projection  of  the  center 
of  the  incident  laser  beam  onto  that  plane  and  pointing  in 
the  direction  cf  the  laser  source. 

The  slopes  of  the  ocean  surface  facet  are  defined  by 
angular  parameters  Alpha  and  Beta.  The  angle  Beta  is  the 
angle  between  the  line  of  steepest  ascent  of  the  facet  and 
the  x-y  plane;  Alpha  is  the  angle  between  the  projection  of 
the  line  cf  steepest  ascent  onto  the  x-y  plane  and  the 
y-axis. 

A  second  cartesian  coordinate  system  at  the  surface  is 
used,  as  explained  in  Ref.  2,  to  simplify  the  computation  of 
the  slope  probability  function.  This  coordinate  system  is 
designated  by  primes  with  the  y '-axis  pointing  towards  the 
wind   source   and  rotated  an  angle  Chi  from  the  y-axis.   The 


z'-axis  coincides  with  the  z-axis.   The  slope  parameters  are 
also  primed  with  their  transformations  given  by: 

Eeta'  =  Beta 
Alpha'  =  Alpha  -  Chi 


Belcw  the  ocean  surface  a  depth  Z  at  the  point  of 
observation,  0',  is  centered  a  translational  transformation 
of  the  surface  coordinate  system.  Angles  Hu  and  Nu  define 
the  laser  beam  refracted  to  0  *  by  any  surface  facet.  The 
first  angle,  Mu,  is  the  angle  between  the  refracted  ray  and 
the  z-axis;  the  second,  Nu,  is  the  angle  between  the  y-axis 
and  the  vertical  plane  of  the  refracted  ray. 


B.   LASER  BEAM 


A  cylindrical  laser  beam  of  radius  a  is  incident  on  the 
ocean  surface  wave  facet  at  an  angle  Psi  measured  from  the 
z-axis.  The  intersection  of  the  beam  and  the  facet  is  a 
flat  elliptical  surface  as  illustrated  in  Fig.  2.  The 
semi-minor  axis  of  the  ellipse  is  equal  to  the  beam  radius 
and  the  semi-major  axis,  b,  is  defined  by: 

b  =  a  *  sec  Psi 


The  analytical  development,  as  in  Refs.  2-4,  begins  at 
the  ocean  surface  facet.  The  laser  source,  aiming, 
atmospheric  propagation,  etc.  are  extraneous  to  this 
analysis.  Pcwer  density  values  at  the  surface  ar€  not 
quoted   because   the   analysis   deals   with  the  ratio  of  the 


power  density  at  the  depth  of  interest,  P  ,   and   the   power 

d 

density  at  the  surface,  P 

ds 


C.   MAXIMUM  ANGLE  OF  INCIDENCE 


The  angle  of  incidence  with  respect  to  the  z-axis,  Psi, 
of  the  laser  team  is  analyzed  from  zero  degrees  to  a  maximum 
of  84  degrees.  At  greater  angles  the  beam  will  not  be  able 
to  strike  all  the  wave  facets  within  a  wavelength  cf  the 
wave  due  to  the  wave  height.  There  is  also  the  possibility, 
as  illustrated  in  Fig.  3,  that  a  ray  incident  at  Psi  >  84 
degrees  could  pass  through  the  wave,  re-enter  the  atmosphere 
and  then  re-enter  the  ocean;  or,  depending  on  the  local 
slopes,  the  ray  could  undergo  a  total  internal  reflection  if 
the  Critical  Angle  were  exceeded.  The  mathematical  model  of 
5wennen[ Ref s.  2-4]  and  this  thesis  do  not  address  this 
situation. 

This  limitation   on   the   maximum   angle   of   incidence, 

Psi    ,  is  derived  from  Fig.  4  and  defined  by: 
max 


-1 
Psi    =  -can   (L/2H) 
max 


The  significant  wave  height,  H,  and  the  average   period,  T, 

of  the  wave  are  given  by  equations  (1)  and  (2)  as  defined  by 

Pierson[Hef.  7].   The  average   wavelength   of   the   wave  is 
given  by  equation  (3)  as  defined  by  Hill[Ref.  8]. 


a  =  2.  14x10    *  V  (1) 


T  =  0.81  *  (2W/g) 


L  =  gT  /2 


(2) 
(3) 


where 

W  =  wind  velocity 

g  =  gravitational  constant 

For   wind  velocities  of  1  to  14  m/sec  the  maximum  angle 

of  incidence  is  84.2  degrees  (±0.4   degrees)   for   all   wind 

velocities;  therefore,  Psi     is  limited  to  84  degrees. 

max 
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III.   POtfER-DENSITY  PROABILITY  INTEGRAL 


A.   GENERAL  IQOATION 


The  power-density  probability  distribution  (P  )   at   the 

d 

point  of  observation  below  the  ocean  surface  (01)  due  to  the 

entire  beam  is  given  by  the  integral  of  dP   over  all  of   the 

d 

contributing   surface   facets.   As  a  function  of  the  angular 

parameters  (tfu,Nu)  the  integral  as  defined  in  Refs.   2-4   is 
given  by: 


P   =  J JF  (Mu,Nu)  dMu  dNu 


(4) 


where 


and 


F(Mu,Nu)  =  P   *T*DTF*cos(WI)*sec  (WE) 

ds 


*[cos  (Mu)  *cos  (Mu  )  +sin  ('An)  *sin  (Mu  )*cos|Nu-Nu  |] 

0  0  0 

2         -1 
*P  (Z'jZ1)  *tan (Beta) *sec  (Beta) *J 
x   y 


P   =  power  density  at  the  surface  facet 
ds 

T  =  Fresnel's  transmittance  coefficient 

DTF  =  diffuse  transmittance  function 

WI  =  angle  of  incidence 
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WE  =  angle  of  refraction 

Mu   =  first  quadrant  angle  the  refracted  center 
0 

ray  of  the  beam  makes  with  the  positive 

z-axis 

Nu   =  positive  angle  the  projection  of  the 

refracted  center  ray  of  the  beam  onto  the 

xy-plane  makes  with  the  y-axis  measured 

from  the  latter  to  the  former 

P(Z*,Z*)  =  time-average  slope  distribution 
x   y 

function 


The  development  of  equation  (4)  is  described  in  great 
detail  in  Refs.  2-4  and  will  not  be  repeated  here;  however, 
for  clarity  the  major  components  of  equation  (4)  will  ne 
discussed. 


B.   FRESNEL'S  TR ANSMITT ANCE  FUNCTION 


Fresnel's  transmittance  coefficient.  (T)  is  defined  as 
the  intensity  ratio  of  the  transmitted  to  the  incident  beam 
energy. 

The  ccllimated  incident  laser  beam  is   composed   of   two 

polarization   components,   a  transverse  magnetic  or  parallel 

polarization  (T   )   and   a   transverse   electric   or   normal 

// 

polarization   (T  )  .    Resolving   the  refracted  beam  into  its 

x 

components  (T  ,T   )  is  a  complex  task  due  to  rotation  of  the 
J.  // 

plane  of  incidence  about  all  three  coordinate  axis  over  the 
range  of  integration.  In  order  to  deal  with  this  complexity 
the  beam  is  assumed  to  be  rendered  diffuse  and  unpolarized 
below  the   ocean   surface   due   to   scattering,   propagation 
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direction   variations   and   phase   changes.   This  assumption 

allows  for  the  averaging  of  T   and  T    to  obtain   the   value 

1      // 

of  I. 


T  =  (T  +T   )/2 
1   // 


The    eguations    for   T      and    T        used    by    S wennen[ Ref s .    2-4] 

1  // 

are: 


T      =    [sin(WI)  *sin  (WR)/sin  (WI+WR)  ] 


(5) 


I         =    [sin  (MI)  *sin  (HE)  /sin  (WI+WR)  *cos  (WI-WR)  ]  (6) 


Eguations  (5)  and  (6)  are  in   disagreement   with   those 
developed  by  Ecrn[5ef.  9]  and  Fowles[Ref.  10],  given  by: 


T   =  [  2*cos(WI)  *sin(WR) /sin  (WI  +  WR)  ] 


(7) 


T    =  [  2*cos(WI)  *sin(WR)  //sin  (WI+WR)  *cos  (WI-WR)  ]     (8) 


Equations  <5)  and  (6)  will  give  erroneously  high  values 
for  the  Fresnel  transmittance  function,  especially  at  higher 
angles  of  incidence.   Because  of  this,  eguations  (7)  and  (8) 
are  used  in  this  modeling. 

The  geometry  of  the  ray  refraction   through   the   ocean 
surface  is  illustrated  in  Fig.  5. 
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DIFFUSE  TRANSMITTANCE  FUNCTION 


The  diffuse  tr ansmittance  function  (DTF)  accounts  for 
beam  attenuation  below  the  ocean  surface  due  to  scattering 
and  absorbticn. 

The  six-ccnstant  DTF  developed  by  Duntley  in  Ref.  1 1  is 
simplified  to  a  function  of  two  constants;  the  backward 
scattering  coefficient  (BS)  and  the  total  absorbtion 
coefficient  (AC) .  This  simplification  is  valid  because  of 
the  randcm  scattering  particle  orientation  in  the  ocean. 
The  DTF  is  defined  by: 


DTF  =  K/[  (AC+BS) *sinh (K*Z) +K*cosh  (K*Z)  ] 

where 


1/2 


/ 

K  =  [  AC*  (AC+2BS)  ] 

Z  =  depth  of  the  point  of  observation 


In  current  terminolgy  this  method  of  accounting  for 
attenuation  is  a  zero-angle  forward  scattering  technique. 

Typical  values  for  BS  and  AC  were  used  from 
measurements  by  Tyler[Ref.  12]  and  Duntley[Ref.  13].  These 
represent  an  average  between  cool  and  warm  ocean  waters  in 
the  blue-green  spectrum  (4800  angstrom). 

BS  =  0.065  percent 

-1 
AC  =  0.044  m 
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D.   SLOPE  PECEABILITY  FUNCTION 


The   time-average   slope  probability  function,  P(Z',Z'), 

x   y 

is  a  statistical  distribution  of  the   ocean   surface   slopes 

developed  by  Cox  and  Hunk[Ref.  5],  The  distribution  derived 
frcm  the  surface  geometry  is  Gaussian,  which  is  then  altered 
by  a  Gram-Charlier  series  in  order  to  account  for  slope 
skewness  and  peakedness  caused  by  the  wind. 


P  (Z*  ,Z»)  =  f  (Xi,Eta,W) 
x   y 


where 


Xi  =  standardized  crosswind  slope  component 
Eta  =  standardized  upwind  slope  component 
w  =  wind  velocity 

The  development  and  the  equation  for  the  slope 
probability  function  are  covered  in  detail  in  Refs.  2,4,  and 
5.  A  clean  ocean  surface  is  assumed  and  the  limits  of 
applicability  placed  on  the  functional  parameters  are 
adhered  tc,  namely: 

| Xi |  <  2.5 

I  Eta j  <  2.5 
W  <  14  m/sec 


E.   JACOEIAN 


The  power-density  probability   integral,   equation   (4) , 
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was   first   developed  over  the  slope  components  (Z1  ,Z%  ) .   It 

x   y 

was  then  transformed  to  a  function  of  the  angular  parameters 

(Alpha, Beta)   and   then  by  the  Jacobian  (J)  to  a  function  of 
(Mu,Nu)  . 


where 


-fl 


dP 


Z'Z' 


r  Tf  (Alpha, Beta)     dAlpha   dBeta 


-1 


f  f  F  (Mu,Nu)     J       dMu    dNu 


J   =    &(Mu,Nu)/^  (Alpha,  Beta) 
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IV. 


METHODS  OF  SOLUTION 


A.   FAR  ZONE 


Swennen[ Rets.  2-4]  observed  that  when  the  depth  (Z)  of 
the  point  of  observation  is  large  compared  to  the  beam 
cross-sectional  radius  at  the  surface  (Z/a>100)  the 
integrand  of  equation  (4)  remains  essentially  constant 
during  the  integration.  The  integral  can  then  be 
approximated  by: 


=     jjF(Mu, 


-1 

Nu) *J         dMu   dNu 


-1 

-   F(Mu,Nu)*J        aHu   4Nu    ;     with    Mu=Mu      &    Nu  =  Nu 

0         0 


The  area  aMu^Nu  in  the  (Mu,Nu)  plane  is  approximated  by 
an  ellipse  of  area  aA,  whose  development  is  covered  in 
detail  in  Ref.  3. 

This  is  the  Far  Zone   approximation   and   equation   (4) 
becomes: 


P   =  p   *T*tan  (Beta) *sec  *cos(WI) 

d    ds 


-1 
*sec  (WR)  *P  (Z»  ,P')  *J   *AA 
x   y 
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The   tasic   Far   Zone   computer   program   developed   in 

FORTRAN   language   and   run  on  the  IBM  360  system  solves  for 

the  ratio  of  the  power  density  at  the  point   of   observation 

to  the  power  density  at  the  surface  facet,  P  /P   .   The  main 

d   ds 

portion   cf   the   program  allows  for   the   entry   of    any 

combination   of   variables   (wind,   depth,  beam  radius,  Psi, 

Chi,   Nu  ,   backward   scattering    coefficient,    absorbtion 

coefficient)   and   also   computes   the  diffuse  transmittance 

function,  Fresnsl* s  transmittance  function  and  the  slope 
probability  function.  Subroutines  common  to  all  variable 
entries  are  used  to  compute  the  angles  (Alpha,  Beta,  HI, 
WS) ,  the  Jaccbian  and   A. 

It  is  possible  to  investigate  the  power-density 
probability  distribution  over  a  wide  range  of  variables  with 
a  minimal  cf  storage  reguirements  and  computation  time  due 
to  the  use  of  an  IBM  System/360  Source  Library  subroutine, 
NLNSYS,  in  solving  for  the  coupled  (Alpha, Beta)  and  (Mu,Nu) 
angle  pairs  from  the  simultaneous  non- linear  eguations: 


cos  (Mu  )  =  [cos  (Psi)  +  K*cos  (Beta)  ]/n  (9) 

0 

cot  (Nu  )  =  cot(Alpha)  -  sin (Psi) /K 
0 

*sin  (Alpha) *sin  (Beta)  (10) 


wnere 

n  =  relative  index  of  refraction  (1.33) 
K  =  K  (Psi, Alpha, Beta, n) 
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B.   NEAfi  ZCNZ 


When  the  depth  of  the  point  of  observation  is  not  large 
compared  to  the  beam  cross-sectional  radius  at  the  surface 
(Z/a<100)  the  Far  Zone  approximation  does  not  hold  and  the 
double  integration  of  equation  (4)  must  be  carried  out. 

Swennen[Ref.  2]  used  Simpson's  method  of  numerical 
integration  to  approximate  equation  (4) .  The  method  was 
subject  to  singularity  points  in  the  integration  limits 
which  increased  the  complexity  and  time  required  for  the 
computation . 

The  doufcle  integration  can  be  more  accurately  and 
relatively  simply  approximated  by  using  the  Gauss  Quadrature 
system  of  solution,  fiefs.  14  and  15.  The  Two-Point 
Gauss-Legendre  Quadrature  method  was  used.  This  method, 
valid  for  polynomials  up  to  degree  three  (equation  (4)  is  of 
degree  two),  consists  of  first  transforming  the  function 
F(Mu,Nu)  into  a  function  F(s,t)  whose  interval  is  -1<s<1  and 
-1<t<1  by  letting: 


Mu  =  [  (Mu  -Mu  )  *s+Mu  +  Mu  ]/2  (11) 

u   1      u    1 

Nu  =  [  (Nu  -Nu  )*t+Nu  +Nu  ]/2  (12) 

u    1       u    1 


and 


dMu  =  [  (Mu  -Mu  ) /2  ]ds  (13) 

u    1 

dNu  =  [  (Nu  -Nu  )  /2  ]dt  (14) 

u    1 
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The  power-density  integral 

P  -   f JP  (Mu,Nu)  dMu  dNu 

is  transformed  into 


il 


P   =  [  (Hu  -Mil  )  (Nu  -Nu  )/4  ]*  f  Tf  (s,t)  ds  dt    (15) 
d        u    1     u    1       J  i 

The      Twc-Point      Gauss-Legendre   Quadrature   approximation 
of    equation    (13)     is   given   by: 

1     1 


P      = 
d 


[  (Mu    -Mu    )   (Nu    -Nu    )/1]*£S»    *F(5    ,t    )        (16) 
u        1  u        1  jso  jso  i  i      J 


where 


W   =  1.0  =  Gauss  weight  factors 
i 


-2 
s  ,t   =  ±  (3)    =  Gauss  roots 
i   J 


The   basic   Near   Zone      computer      program      developed      uses 

this        Gauss-Legendre        Quadrature      method      in      solving      the 

power-density   integral   coupled    with   the   same   subroutines      as 

the      Far      Zone      program.         The      main      program      computes      the 

integration    limits    (Mu   ,Mu   )     and    (Nu    , Nu   )     after   determining 

u    1         u   1 

whether   the   surface  projection  of  the  point  of  observation 

is  inside  or  outside  the  beam's   horizontal   cross-secticnal 

area.    The   equations   for   the  integration  limits  are  well 

defined  in  Refs.  2-4. 

The  transformations  described  by  equations  (11)- (15) 
are  than  performed  and  the  power-density  is  solved  for  using 
equation   (16).    As   in  the  Far  Zone  program,  the  Near  Zone 
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program  can  take  any  combination  of  variables  as  entries 
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7.       RESULTS 


Results      for   oblique   incidence   angles    (0<Psi<84    degrees) 

were   obtained    for    vertical    plane   cuts   at      Nu        =0      degrees 

0 

over   the  entire  range  of  variables;  however,  valid  data  for 

vertical  planar  sections  at  Nu   #0   degrees   could   not   be 

0 

obtained   due   to   values   of  the  standardized  crosswind  and 

upwind   slope   components   being   out   of    the    range   of 

applicability   for   the   slope   probability   function.   This 

problem   is   believed   to   be   inherent   in   the   method   of 

mathematical   analysis   of   the  ocean  surface  used.   Results 

for  normal  incidence  (Psi  =  0  degrees)  over  the  entire  range 

of    variables,    including    Nu   planar   variations,   were 

0 

obtained. 


The  validity  check  carried  out  on  the  slope  probability 
function,  as  discussed  in  Section  III.D,  clearly  identified 
invalid  data.  Another  source  of  invalid  data  occurred  at 
higher  incidence  angles  (Psi)  where  the  simultaneous 
solution  of  ncn-linear  equations  (9)  and  (10)  resulted  in  an 
angle  of  incidence  (WI)  greater  than  90  degrees  while  the 
slope  probability  function  still  indicated  valid.  This 
output  was  also  easily  identified  and  eliminated. 

Computations  with  the  Near  Zone  solution  method  used 
more  computer  time  and  storage  than  did  the  Far  Zone 
solution  method.  For  a  computer  solution  using  one  value 
for  each  variable  entry,  the  Near  Zone  used  36  per  cent  more 
computer   time  and  8  per  cent  more  computer  storage  than  did 
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an  equivalent  Far  Zone  solution.  The  higher  time 
requirement  for  the  Near  Zone  solution  was  due  mainly  to  a 
greater  compiling  time  requirement.  The  time  difference 
between  the  two  solution  methods  decreased  as  the  number  of 
variables  that  were  incremented  during  one  computer  run 
increased. 

A  sample   computer  run   with   output   is   presented   in 
Appendix  E. 


A.   PRESENTATION  OF  RESULTS 


Results  are  presented  graphically  to  simplify 
comparisons  cf  the  wide  range  of  variables.  Vertical-planar 
sections  cf  the  power  density  probability  distribution  below 
the  ocean  surface  cut  through  the  center  of  the  beam 
cross-section  at  the  horizontal  ocean  surface   at   an   angle 

Nu    are   used   for   the   majority  of  the  results  presented. 
0 

This   presentation   is  in   rectangular   coordinates   with 

normalized  power  density  in  decibels  on  the  ordinate  and  the 

angle  Mu   in  degrees  on  the  abscissa.   The  normalized   power 
0 

2 
density     is    equal   to    10*log   [ (P  *Z  ) /P   ].     This 

10    d      ds 

2 
normalization   with   respect   to   P  /Z    should   produce   a 

ds 

distribution   function   that   is  independent  of  depth  in  the 

Far  Zone  regime,  according  to  Swennen[Ref.  2-4],  because  the 
multiple  slopes  seen  in  a  beam  cross-section  at  the  surface 
tend  to  render  the  radiation  diffuse.  The  power  density  of 
diffuse  radiation  at  large  depths  decreases  as  the  square  of 
the  depth.   When  the  effects  of   scattering   and   absorption 
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are  included  this  independence  does  not  occur,  as  will  be 
discussed  later;  however,  the  normalization  of  the  power 
density  is  used  in  all  regimes  for  continuity  and  for 
comparisons  with  Refs.  2-4. 


E.   EFFECTS  CF  WIND  VELOCITY 


As  wind  velocity  increases,  the  power  density 
distribution  spreads  out  and  the  maximum  power  density 
decreases  slightly.  This  effect  is  presented  in  Figs.  6-11 
for  various  angles  of  beam  incidence  with  respect  to  the 
vertical  (Psi  =  0,30,60  and  84  degrees).  Figures  10  and  11 
also  show  the  close  correlation  between  results  obtained 
from  the  Near  Zone  and  the  Far  Zone  solution  methods  for  Psi 
=  0  and  40  degrees,  respectively. 

The  spreading  effect  on  the  power  density  distribution 
is  due  to  larger  refraction  angles  (WR)  caused  by  higher 
wave  facet  slopes  that  occur  at  higher  wind  velocities.  The 
maximum  power  density  decrease  is  caused  by  the  spreading. 

For  an  optical  communication  or  detection  system  an 
increase  in  wind  velocity  would  mean  that  the  maximum 
intensity  cf  the  beam  that  could  be  focused  to  a  point 
beneath  the  ocean  surface  would  be  reduced,  but  the  width  of 
receivable  signal  radiation  or  the  search  width  would  be 
increased. 


C.   EFFECTS  OF  INCIDENCE  ANGLE 


As   the   angle  of  incidence  with  respect  to  the  vertical 
(Psi)  is  increased  the  maximum  power  density   decreases   and 
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the  distribution  shifts  away  from  the  vertical.  This  effect 
can  be  seen  in  Figs.  6-11  and  more  clearly  in  Figs.  12  and 
13.  Figure  12  is  a  Far  Zone  solution  at  a  depth  of  50 
meters  for  Psi  =  0,30,60  and  84  degrees  while  Fig.  13  is  a 
Near  Zone  solution  at  a  depth  of  10  meters  for  the  same 
values  of  Psi. 

The   decrease   in   maximum   power  density  is  a  result  of 

greater  reflectance  at  the  surface  facet  as   the  angles  of 

incidence  (WI)  increase  with  increases  in  Psi.   The  shift  in 

the  center  of  the  distributions  away  from   the   vertical  is 

caused  by  the  higher  angles  of  refraction  (HE)  that  occur  as 
WI  increases. 

For  an  optical  communication  or  detection  system  an 
increase  in  Psi  would  be  expected  to  lower  the  maximum  beam 
intensity  below  the  ocean  surface  and  to  shift  the  point  of 
maximum  intensity  away  from  the  point  of  water  entry. 


D.   EFFECTS  CF  HIND  DIRECTION 


The  effects  of  wind  direction  are  small  compared  to  the 
wind  velocity  effects.  When  the  wind  direction  (Chi)  is 
perpendicular  to  the  beam  (Chi  =  90  degrees)  there  is 
slightly  less  spreading  of  the  power  density  distribution 
compared  to  Chi  =  0  or  180  degrees.  The  maximum  power 
density  is  essentially  unaffected  by  variations  in  Chi. 

Figure  14  shows  the  effect  of  variations  in  Chi  for 
normal  incidence  and  Fig.  15  for  oblique  incidence  (Psi  =  40 
degrees) . 

Wind  direction  appears  to  be  of  little  concern  to  an 
optical   communication  or  detection  system  as  related  to  the 
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power  density  distribution  below  the  ocean  surface 


E.   EFFECTS  CF  SPOT  SIZE 


Figures  16  and  17  show  the  effect  of  increasing  the 
incident  beam  spot  radius  from  0.1  to  0.5  meters  for  Esi  = 
0  and  40  degrees,  respectively.  This  is  done  with  the 
assumption  that  the  energy  per  unit  area  remains  constant; 
therefore,  the  total  energy  must  increase.  The  power 
density  spreads  out  and  increases  in  intensity  uniformly 
over  the  distribution  as  the  spot  size  is  increased. 

Since  the  equations  used  in  this   solution   method   deal 

with   the   ratio   of   P  /P   ,  these  results  simply  mean  that 

d   ds 

more  power  in  a  larger  beam  will  produce  a  more  intense   and 

more  widely  distributed  pattern  below  the  ocean  surface. 


F.   EFFECTS  OP  DEPTH 


The  decrease  of  power  density  with  increases  in  depth  is 
illustrated  in  Figs.  18  and  19  f or  Psi  equal  to  0  and  40 
degrees,  respectively.  The  decrease  is  linear  for  the 
normalized  power  density  as  can  be  seen  in  Fig.   20. 

Scattering  and  absorption  cause  this  decrease  in  tne 
power  density.  As  the  depth  increases  the  diffuse 
transmittance  function  (DTF)  rapidly  becomes  an  inverse 
function  of  the  sum  of  the  hyperbolic  sine  of  the  depth  and 
the  hyperbolic  cosine  of  the  depth, 
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^1 

DTF  =  f[  (sinhZ+coshZ)   ] 


which  means  the-  DTF  decreases  rapidly  with  increases  in 
depth  and  approaches  zero  in  the  limit.  This  agrees  with 
the  theory  cf  Preisendorf er[ Ref .  16]. 

Swennen's  prediction  of  a  constant  normalized  power 
density  distribution  in  the  Far  Zone  regime  due  to  the 
diffusing  effects  of  the  ocean  in  Refs.  2-4  holds  only  when 
the  effects  of  scattering  and  absorption  are  ignored.  One 
computer  run  was  made  with  the  DTF  set  equal  to  one  to 
illustrate  this  effect  in  Fig.  20. 

Figures  21  and  22  illustrate  the  decrease  in  power 
density  with  depth,  non-dimensionalized  with  respect  to  spot 
size,  Z/a,  again  for  Psi  equal  to  0  and  40  degrees, 
respectively. 

An  optical  communication  or  detection  system  will  be 
limited  greatly  by  the  power  density  decrease  associated 
with  depth. 


G.   MAXIMUM  POWER  DENSITY 


The  maximum  obtainable  power  density  available  over  the 
spectrum  of  incidence  angles  (Psi  =0-84  degrees)  for  a 
given  set  of  conditions  is  presented  in  Figs.  23  and  24. 
The  maximum  power  density  obtainable  does  not  occur  for 
normal  incidence,  but  for  a  small  value  of  Psi,  because  of 
the  fact  that  as  the  wind  increases  the  most  probable  slope 
is  not  zero  but  a  small  angle.  The  maximum  obtainable  power 
density  then  drops  continuously  with  increases  in  Psi  after 
the  peak  that  occurs  around  Psi  =10-15  degrees. 
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Figure  23  illustrates  the  decrease  in  the  maximum 
obtainable  power  density  for  all  values  of  Psi  as  the  wind 
velocity  increases.  This  is  due  to  the  spreading  effect  on 
the  power  density  distribution  caused  by  the  wind  that  was 
described  previously. 

Figure  24  illustrates  the  independence  of  the  maximum 
power  density  on  the  wind  direction  (Chi)  for  all  values  of 
Psi.  Chi  equal  to  270  degrees  is  not  plotted  but  is 
equivalent  to  the  condition  at  90  degrees  because  of  the 
symmetry  of  the  slope  probability  function  about  the  wind 
direction. 

The  angle  with  respect  to  the  vertical  from  the  point  of 

observation  helow  the   ocean   surface   (Mu  )   at   which   the 

0 

maximum   power   density   occurs   for   each   value   of  Psi  is 

plotted  in  Fig.    25.    The   plot   is   independent   of   wind 

velocity   and   direction.   Examining  Figs.  16  -  19  it  can  be 

seen  that  Fig.  25  is  also   independent   of   depth   and   spot 

size.    Mu    for   maximum   power   density   increases   nearly 
0 

linearly  frcm  approximately  1  to  56  degrees  as  Psi  goes  from 

0   to   84   degrees.   The  increase  in  Mu   is  due  to  the  power 

0 

density  distribution  shifting  away   from   the   vertical   for 

increasing    valules   of    Psi   as   discussed   earlier   and 
illustrated  in  Figs.  12  and  13. 

The   obvious  consideration  with  an  optical  communication 

or  detection  system  is  the  decrease  in  maximum  power  density 

that    accompanies    higher    incidence    angles    and   wind 

velocities.   The  strong  independence  of   the   angle   Mu   at 

0 

which   the   maximum  power  density  occurs  with  respect  to  the 
wind  velocity,  wind  direction,  depth  and  spot  size  indicates 
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that  a  satellite  based  laser  navigation  system  for  submarine 
usage  may  be  a  possibility. 


29 


71.   CONCLUSIONS 


The  time-averaged  power-density  distribution  below  the 
ocean  surface  due  to  incident  laser  radiation  has  been 
examined  through  computer  simulations  in  both  the  Near  Zone 
and  the  Far  Zcne  regimes.  The  effects  of  altering  incidence 
angles,  wind  velocities,  wind  directions,  beam  spot  sizes 
and  depths  of  penetration  have  been  presented. 

The  power-density  distribution  was  found  to  be  highly 
dependent  on  the  angle  of  beam  incidence,  the  wind  velocity 
and  the  attenuation  associated  with  depth  of  penetration  of 
the  laser  beam.  The  fact  that  the  location  of  the  maximum 
power  density  was  found  to  be  fundamentally  dependent  on 
only  the  angle  of  beam  incidence  indicates  a  possibility  of 
direction  finding  capabilities  by  a  submerged  receiver. 

Improvements  in  computer  systems  and  analysis  techniques 
have  allowed  an  examination  of  a  much  broader  range  of 
variables  than  were  possible  heretofore.  Angles  of  beam 
incidence  up  to  84  degrees  were  considered  within  the  range 
of  validity  fcr  restrictions  placed  on  the  simulation. 
Above  84  degrees  incidence  angle  the  possible  interactions 
at  the  air/sea  interface  make  an  analysis  of  the 
power-density  distribution  extremely  complex,  and 
experimentation  will  be  necessary  to  give  indications  of  the 
feasibility  of  laser  operations  in  this  area. 

This  simulation  could  be  extended  to  an  analysis  of  much 
larger  spot  sizes  that  would  be  expected  from  a  satellite 
based  transmitter.  Prettyman  and  Cermak[Ref.  17]  suggested 
that   the   effects   of  larger  spots  could  be  approximated  by 
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assuming  €ach  wave  facet  independent  of  its  neighbors  and 
then  integrating  over  the  larger  spot  area,  since  the  small 
spot  represents  the  outer  limit  of  the  effects  of  the  ocean 
surface  on  a  laser  beam  of  greater  spot  size.  There  would 
have  to  be  a  large  number  of  surface  facets  contained  in  the 
spot  to  insure  facet  independence  and  accurate 
approximations,  hence  a  strong  dependence  on  the  wavelength 
of  the  wave  associated  with  wind  velocity. 

An  extension  of  the  simulation  to  pulsed  laser  radiation 
is  also  warrented,  since  this  mode  of  operation  may  be 
reguired  for  air-to-subsurface  laser  systems.  The  effects 
of  pulsed  radiation  occur  nearly  instantaneously  as  compared 
to  the  time-averaged  analysis  performed  in  this   simulation. 

A  further  refinement  of  the  simulation  would  involve  a 
more  complex  modeling  of  the  beam  attenuation  due  to 
scattering  and  absorption  in  the  ocean  through  the  use  of  a 
diffusion  or  a  multiple  forward  scattering  method,  both 
described  by  Bucher[Sef.  18]  and  Preisendorf er[ Ref .  16]. 


31 


APPENDIX  A. 
FIGURES 
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Figure  1:  Coordinate  systens 
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Figure   2:      Intersection  of   laser  beam 
with   ocean  wave    facet 
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Figure  3:      Possible   beam  paths    at 
high   incidence   angles 
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Figure  4:   Maximum  angle  of  incidence 
wave  geometry 
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Figure   5:      Ray  refraction  geometry 
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APPENDIX  B 


COMPUTER  PROGRAM  USAGE 


To  use  the  computer  programs  listed  following  this 
Appendix,  first  compute  the  ratio  of  the  depth  of  the  point 
of  observation  below  the  ocean  surface  to  the  beam  radius  at 
the  surface  facet,  Z/a.  Use  the  Near  Zone  program  for 
Z/a<100  and  the  Far  Zone  program  for  Z/a>100.  The  Near  Zone 
program  may  be  used  for  any  value  of  Z/a,  but  a  slight 
computation  time  and  computer  storage  increase  will   result. 

The  only  other  basic  requirement  is  the  selection  of 
eight  input  variables.  The  variables  are  listed  below  along 
with  the  line  from  the  respective  program  in  which  they 
appear. 

NEAR   FAR 
VARIABLE  ZONE   ZONE 

Depth,  Z  (m)  470    440 

Beam  radius,  RAD  (m)  480    450 

Angle  Nu  ,  NUORUN  (deg)      .  490    480 

0 

Backward  scattering  coefficient,  BS        560    530 

-1 

Absorption  coefficient,  AC  (m   )  570    540 

Beam  incidence  angle,  PSIRUN  (deg)         630    590 
Wind  direction,  CHI  (rad)  680    640 

Wind  velocity,  W  (m/sec)  710    670 

Z,   RAD,   NUORUN,   BS   and  AC  are  direct  program  inputs 
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made  on  the  indicated  lines.  PSI,  CHI  and  W  are  selected 
from  DATA  statements  430  and  440  in  the  Near  Zone  program 
and  420  and  430  in  the  Far  Zone  program.  The  variable 
locations  in  DATA  statements  were  used  for  multiple 
variations  of  the  parameters  in  a  slightly  different  form  of 
the  program  than  presented  which  is  easily  set  up  with  DO 
LOOPS.  The  entry  variables  may  also  be  varied  in  this  same 
manner. 

Both  programs  must  be  combined  with  the  common  package 
of  subroutines  that  are  listed  following  the  basic  program 
listings. 

The  output  from  both  listed  programs  appear  on  the  next 
two  pages. 
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